Introducción

Descargar los datos:

Para descargar los datos, es importante conocer en que extención se encuentran guardados y así utilizar el recurso correcto. El siguiente es el utilizado par datos guardados con excel.

readxl: Necesario para leer archivos excel

library(readxl)

Descarga de base de datos:

Cuando se descargue la base de datos, es importante tener encuenta que:

1. Se tiene información histórica suficiente dependiendo de la periodicidad definida al inicio del estudio. Esto quiere decir que si se tiene información mensual debería ser por varios años, si se determina que es un estudio de algunos meses, debe tenerse por varios días.

2. Esta información debe ser continua en el tiempo, es decir que cubra todo el periodo escogido y que la separación temporal sea constante.

3. Que si falta información sea muy poca, preferiblemente no mayor al 20%, lo que permitiría calcular los datos faltantes para que no comprometa de forma importante la información que se pueda extraer de los datos.

Ejemplo: Estos datos corresponden a veinte años (enero de 1979 a diciembre de 1998) de registros de muertes totales en Colombia, que se encuentran en la pagina del DANE (https://www.datos.gov.co/widgets/kk5w-ugzm), es decir 240 datos. Como se puede ver, la base de datos tiene tres columnas, la primera de ella es un número consecutivo desde la fecha más lejana en el tiempo, hasta la fecha más reciente, la segun hace referencia a los meses y la tercera columna son los conteos de muertes por mes en cada año.

ColombiaT <- read_excel("ColombiaT.xlsx")

## New names:
## • `` -> `...1`

head(ColombiaT)

## # A tibble: 6 × 3
##   ...1  MES   conteo
##   <chr> <chr>  <dbl>
## 1 1     Ene     9582
## 2 2     Feb     8467
## 3 3     Mar     9202
## 4 4     Abr     9197
## 5 5     May     9484
## 6 6     Jun     9596

Revisión de base de datos

Esto sirve para identificar las variables de la base de datos y que tipo de variables son.

summary(ColombiaT)

##      ...1               MES                conteo     
##  Length:240         Length:240         Min.   : 6980  
##  Class :character   Class :character   1st Qu.:11641  
##  Mode  :character   Mode  :character   Median :12842  
##                                        Mean   :12760  
##                                        3rd Qu.:14051  
##                                        Max.   :18507

Datos faltantes

Librerias necesarias:

mice: Este paquete lo ayuda a imputar valores faltantes con valores de datos plausibles. Estos valores plausibles se extraen de una distribución diseñada específicamente para cada punto de datos faltante.La creación de múltiples imputaciones para que sea seleccionada la más adecuada.

VIM: Herramientas para la visualización de valores perdidos y/o imputados, que pueden utilizarse para explorar los datos y la estructura de los valores perdidos y/o imputados. Dependiendo de esta estructura de los valores perdidos, los métodos correspondientes pueden ayudar a identificar el mecanismo que genera los valores perdidos y permite explorar los datos, incluidos los valores perdidos. Además, la calidad de la imputación se puede explorar visualmente utilizando varios métodos de gráficos univariados, bivariados, múltiples y multivariados..

missForest: La función ‘missForest’ en este paquete se usa para imputar valores faltantes, particularmente en el caso de datos de tipo mixto. Utiliza un bosque aleatorio entrenado en los valores observados de una matriz de datos para predecir los valores que faltan. Se puede utilizar para imputar datos continuos y/o categóricos, incluidas interacciones complejas y relaciones no lineales. Produce una estimación del error de imputación out-of-bag (OOB) sin la necesidad de un conjunto de pruebas o una validación cruzada elaborada. Se puede ejecutar en paralelo para ahorrar tiempo de cálculo..

library(mice)

## 
## Attaching package: 'mice'

## The following object is masked from 'package:stats':
## 
##     filter

## The following objects are masked from 'package:base':
## 
##     cbind, rbind

library(VIM)

## Loading required package: colorspace

## Loading required package: grid

## VIM is ready to use.

## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues

## 
## Attaching package: 'VIM'

## The following object is masked from 'package:datasets':
## 
##     sleep

library(missForest)

## 
## Attaching package: 'missForest'

## The following object is masked from 'package:VIM':
## 
##     nrmse

En el analisis de series de tiempo, un requerimiento es que no existan datos faltantes. En esta base de datos noa hay datos faltantes, pero para explicar como proceder en este caso, generaremos algunos en la base de datos y procederemos a calcularlos, lo que se conoce como imputación de datos.

Ejemplo de imputación

Como ejemplo para mostrar como proceder cuando se encuentran valores ausentes, se tomarán los datos y se eliminarán algunos valores, para luego remplazarlos por unos generados automáticamente con algunas condiciones. Se generá una base de datos llamada ColombiaNA con 5% de valores ausentes, que se seleccionaron de forma aleatoria, y corresponden a 36, como se puede ver en el summary de los datos

Antes de eliminar algunos valores, adicionaremos dos variables que tiene una distribución normal, para tener tres varibles numéricas.

set.seed(123)
Colombiamiss<-ColombiaT
Colombiamiss$Variable1<-rnorm(240)
Colombiamiss$Variable2<-rnorm(240)

Aquí lo que se hace es eliminar 5% en cada variable de forma aleatoria, para tener NAs y mostrar el procedimiento para la imputación. Se tiene un set.seed, para asegurar que siempre que se eliminen datos no varien los resultados. Esto crea resultados reproducibles completamente. Al final se produce un resultado que es el resumen donde se pueden ver las características de las variables y la cantidad de NA en cada una de ellas.

set.seed(123)
ColombiaNA<-prodNA(Colombiamiss[,-c(1,2)],noNA = 0.05)
summary(ColombiaNA)

##      conteo        Variable1           Variable2       
##  Min.   : 8424   Min.   :-2.309169   Min.   :-2.66092  
##  1st Qu.:11636   1st Qu.:-0.637506   1st Qu.:-0.55858  
##  Median :12824   Median :-0.086844   Median : 0.07455  
##  Mean   :12760   Mean   :-0.004563   Mean   : 0.06539  
##  3rd Qu.:14045   3rd Qu.: 0.579208   3rd Qu.: 0.78722  
##  Max.   :18507   Max.   : 3.241040   Max.   : 2.57146  
##  NA's   :13      NA's   :10          NA's   :13

Identificación de valores ausentes Podemos usar md.pattern() para tener una idea de los patrones de datos faltantes. Esto nos muestra que en 205 filas no hay datos faltantes en ninguna variable (la fila superior azul); la siguiente fila dice que faltaban 12 individuos solo en la variable tercera, mientras que en las otras variables los datos estan completos; la fila tres muestra que faltan 12 datos en la segunda variable únicamente. Así es la forma de leer la figura que a continuación aparece. Finalmente, se encuentra que faltan 36 datos en toda la base considerada, que es el valor que aparece en la esquina inferiro derecha.

md.pattern(ColombiaNA)

## Error : The fig.showtext code chunk option must be TRUE

##     Variable1 conteo Variable2   
## 205         1      1         1  0
## 12          1      1         0  1
## 12          1      0         1  1
## 1           1      0         0  2
## 10          0      1         1  1
##            10     13        13 36

El gráfico de abajo nos ayuda a entender nuevamente los valores faltantes de forma diferente, pero puede ser un recurso importante para ener datos en procentajes, utilizando el recurso siguiente: Colplot$percent: aquí encontrará valores en porcentaje de las diferentes conbinaciones de variables, por ejemplo el 85.4% de los valores no tienen datos faltantes en ninguna variable, el 5% solo la segunda variable tiene datos faltantes, el 4% solo la primera variable tiene datos faltantes, el 5% la tercera variable solo tiene datos faltantes y por últimoel 0.4% presenta datos faltantes en las dos primera variables, es decir que comparten daos faltantes esas dos variables.

Colplot<-aggr(ColombiaNA, col=c("blue","red"),
              numbers=TRUE, sortVars=TRUE,
              labels=names(ColombiaNA), cex.axis=.7,
              gap=3, ylab=c("Datos ausentes","Parametros"))

## Error : The fig.showtext code chunk option must be TRUE

## Error : The fig.showtext code chunk option must be TRUE

## Error : The fig.showtext code chunk option must be TRUE

## Error : The fig.showtext code chunk option must be TRUE

## 
##  Variables sorted by number of missings: 
##   Variable      Count
##     conteo 0.05416667
##  Variable2 0.05416667
##  Variable1 0.04166667

El siguiente recurso es una matriz donde se identifican la cantidad de valores faltantes Los cuatro componentes en el valor de salida tienen la siguiente interpretación:

Patrón rr: Se observan ambas variables donde no hay valores faltantes Patrón rm: Se observa la primera variable y falta datos de la segunda variable; Patrón mr: Faltan datos de la primera variable, mientras que en la segunda variable estan completos Patrón mm: Coindiden las dos variables donde no hay datos en ninguna de ellas.

md.pairs(ColombiaNA)

## $rr
##           conteo Variable1 Variable2
## conteo       227       217       215
## Variable1    217       230       217
## Variable2    215       217       227
## 
## $rm
##           conteo Variable1 Variable2
## conteo         0        10        12
## Variable1     13         0        13
## Variable2     12        10         0
## 
## $mr
##           conteo Variable1 Variable2
## conteo         0        13        12
## Variable1     10         0        10
## Variable2     12        13         0
## 
## $mm
##           conteo Variable1 Variable2
## conteo        13         0         1
## Variable1      0        10         0
## Variable2      1         0        13

Esta gráfica pretende mostrar la relación de dos variables, y los datos faltantes quí estamos limitados a graficar 2 variables a la vez, pero, sin embargo, podemos recopilar algunas ideas interesantes. El gráfico de caja roja de la izquierda muestra la distribución de la variable 1 donde presenta los datos faltantes de la variable conteo, mientras que el gráfico de caja azul muestra la distribución de los puntos de datos restantes. Lo mismo ocurre con los diagramas de caja roja de la variable conteo, donde se presentan los valores faltantes de la variable 1 y la caja azul donde los datos estan completos. Si nuestra suposición es correcta, entonces esperamos que los diagramas de caja roja y azul sean muy similares.

marginplot(ColombiaNA[,c("conteo","Variable1")])

## Error : The fig.showtext code chunk option must be TRUE

imputar<-mice(ColombiaNA,m=3, seed = 123)

## 
##  iter imp variable
##   1   1  conteo  Variable1  Variable2
##   1   2  conteo  Variable1  Variable2
##   1   3  conteo  Variable1  Variable2
##   2   1  conteo  Variable1  Variable2
##   2   2  conteo  Variable1  Variable2
##   2   3  conteo  Variable1  Variable2
##   3   1  conteo  Variable1  Variable2
##   3   2  conteo  Variable1  Variable2
##   3   3  conteo  Variable1  Variable2
##   4   1  conteo  Variable1  Variable2
##   4   2  conteo  Variable1  Variable2
##   4   3  conteo  Variable1  Variable2
##   5   1  conteo  Variable1  Variable2
##   5   2  conteo  Variable1  Variable2
##   5   3  conteo  Variable1  Variable2

pmm: Predictive mean matching

print(imputar)

## Class: mids
## Number of multiple imputations:  3 
## Imputation methods:
##    conteo Variable1 Variable2 
##     "pmm"     "pmm"     "pmm" 
## PredictorMatrix:
##           conteo Variable1 Variable2
## conteo         0         1         1
## Variable1      1         0         1
## Variable2      1         1         0

imputar$imp$conteo

##         1     2     3
## 14  11228 11801 11515
## 23  13079 15236 13459
## 26  13905 10663 14019
## 91  11525 13683 12764
## 118 13849 14321 12754
## 135 13079 11515 11133
## 143 12366 14600 15042
## 166 12366  9049 13229
## 179 11525 11394 12168
## 195 14277 12875 11306
## 211 12322 12800 13643
## 224 13918  9049 11554
## 229 14235 15406 14571

Datos completos

DatosF<-complete(imputar,1)

Datos observados vs imputados Se puede obtener otra toma visual útil de las distribuciones usando la función stripplot() que muestra las distribuciones de las variables como puntos individuales.

stripplot(imputar,pch=20,cex=1.2)

## Error : The fig.showtext code chunk option must be TRUE

Comparando la distribución de los datos originales y los imputados. En este caso se utiliza el diagrama de dispersión del conteo con lad sos varialbes restantes 1 y 2. Se puede observar que los valores imputados, aquí en color verde, estan muy cercanos a los valores originales en color amarillo

xyplot(imputar,conteo ~ Variable2+Variable1,pch=18,cex=1)

## Error : The fig.showtext code chunk option must be TRUE

Lo que nos gustaría ver es que la forma de los puntos verde (imputados) coincida con los amarillos (observados). La forma coincidente nos dice que los valores imputados son de hecho “valores plausibles”.

Otro gráfico útil es el gráfico de densidad: La densidad de los datos imputados para cada conjunto de datos imputados se muestra en magenta, mientras que la densidad de los datos observados se muestra en azul. Nuevamente, bajo nuestros supuestos anteriores, esperamos que las distribuciones sean similares.

densityplot(imputar)

## Error : The fig.showtext code chunk option must be TRUE

Para realizar un análisis de series de tiempo es indispensable tener claros varios aspectos tales como:

Objetivo(s) que se buscan
Identificar los datos que se requieren para cumplir ese(os) objetivo(s)
Determinar la periodicidad que se requiere (años, meses, días, horas, etc.)
Tener claro el nivel de cobertura que se va a utilizar (municipal, regional, departamental, nacional, etc.)
No tener datos perdidos(NAs) y si los tiene debe ser menor al 20%, ya que con esta copndición se pueden calcular los perdidos por varios metodos que se conocen como imputación de datos.

Librerias necesarias:

En R se encuentram una gran cantidad de recursos y librerias que permiten trabajar series de tiempo. Los siguientes son los recomendados por una gran cantidad de personas que trabajan en estos temas. Según la Red Integral de Archivos en R (CRAN: por sus siglas en ingles: Comprehensive R Archive Network), los paquetes presentados aquí son: modeltime: Es una extensión del ecosistema Tidymodels para el modelado de series temporales. Los modelos incluyen ARIMA, Suavizado exponencial y modelos de series de tiempo adicionales de los paquetes ‘forecast’ y ‘prophet’.

tidyverse: Es un conjunto de paquetes que permite la organización de bases de datos, el manejo y la graficación de resultados. Los paquetes incluidos son: dplyr,ggplot2, forecast, tibble, readr, stringr, tidyr y purrr.

forecast: Métodos y herramientas para mostrar y analizar pronósticos de series temporales univariadas, incluido el suavizado exponencial a través de modelos de espacio de estado y el modelado ARIMA automático.

tseries: Permite plotear y análizar series de tiempo de forma sencilla.

gtable: Herramientas para facilitar el trabajo con tablas.

TSstudio: Proporciona un conjunto de herramientas para el análisis descriptivo y predictivo de datos de series temporales. Eso incluye funciones para la visualización interactiva de objetos de series temporales y también funciones de utilidad para la automatización de pronósticos de series temporales.

Las figuras generadas por este medio, presentarán una serie de herramientas en la parte superior: La primera de ellas es la imagen de una camara fotográfica que permite descargar la imagen en .png, a continuación encontrará una lupa para realizar zoom a la figura en un sector específico, luego una cruz, que es una herramienta si es necesario navegar por entre la figura debido a que el pantallazo es muy pequeño. A continuación, aparecen un más (+) y un menos (-), que tambien es un tipo de zoom pero sin especificar la region que se quiere revisar. Luego aparece una X en un recuadro que permite volver al estado inicial de la imagen eliminando los zoom previos.

library(modeltime)
library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.1
## ✔ readr   2.1.2     ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks mice::filter(), stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(tseries)
library(gtable)
library(TSstudio)
library(bookdown)

Conversión a serie de tiempo:

Para poder hacer un análisis de series de tiempo debe estar en forma de matriz donde la primera columna son los años y las siguientes columnas son los meses como aparece a continuación. Se recomienda generar las fechas a partir de un script de R, para no tener problemas de reconocimiento por parte del sistema. El script presentado abajo, inicia teniendo en cuenta la tercera columna (conteo), luego se conbierte los datos en series de tiempo con la función ts y se le especifica el rango del tiempo considerado, especificando donde inicia 1979,1(1979, enero) y donde termina 1998-12 (1998, diciembre), con una frecuencia de 12, ya que se tiene información de los doce meses del año. Cuando se pide mostrar la serie de tiempo (Datos.ts), aparece la estructura. Al utilizar la función glimpse(Datos.ts), como resultado aparece la confirmación que los datos ya tienen un formato de serie de tiempo y las funciones start(Datos.ts);end(Datos.ts), nos mostrará la fecha de inicio y la fecha de terminación de la serie de tiempo.

Datos.ts<-ColombiaT$conteo %>% 
  ts(start=c(1979,1), end=c(1998,12),frequency=12)
Datos.ts

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 1979  9582  8467  9202  9197  9484  9596  9560  9049  8424  8988  8715  8899
## 1980  8883  6980 10443  9841 10142  9919 11473 11353 10847 11515 11067 11957
## 1981 12244 10509 11503 11042 12112 11808 12145 11765 11096 11314 11054 11555
## 1982 12008 10478 11160 11339 11394 11478 11804 11626 10922 11849 11825 11909
## 1983 12764 10663 12120 11378 12071 11435 12355 11772 11239 11554 11139 11967
## 1984 12169 11075 11591 11021 11215 11142 11466 11245 11135 11525 11377 12366
## 1985 18507 13849 11759 10896 12100 11801 12387 12236 12168 12345 12757 13252
## 1986 15773 11133 12002 11663 11645 12084 12458 11967 11669 11909 11962 12202
## 1987 12322 11228 12101 11719 12754 12437 13492 12261 12043 12006 14507 15236
## 1988 11290 12861 12916 11306 13079 13035 13229 11628 12382 12723 14986 13823
## 1989 14120 12056 12875 13031 12610 12684 13348 12526 12375 12824 12917 13459
## 1990 13786 12098 13767 13183 13234 12949 13406 13089 12467 12911 12573 11368
## 1991 14277 12381 13112 13399 14101 13578 14121 13847 13781 13635 13339 14281
## 1992 14042 12681 13774 13131 13954 13683 14321 14235 14029 14346 14034 15644
## 1993 15041 13143 13973 13278 14534 14019 14465 14266 13692 13992 13577 14833
## 1994 15038 12597 13544 13849 14600 14447 14441 14135 13643 14060 13905 14484
## 1995 14629 12869 13754 13043 14837 14853 15406 14571 13918 13758 13355 15093
## 1996 14905 13279 14306 14048 14419 14339 14577 14386 14509 15233 14680 15001
## 1997 15042 12800 14204 13892 14699 14403 14706 14430 13536 13604 14234 15390
## 1998 15085 13311 14350 13531 14586 14264 14610 14959 14906 15124 15105 15703

glimpse(Datos.ts)

##  Time-Series [1:240] from 1979 to 1999: 9582 8467 9202 9197 9484 ...

start(Datos.ts);end(Datos.ts)

## [1] 1979    1

## [1] 1998   12

Análisis exploratorio :

Componente estacionario

El recurso inicial visual de la presentación de los datos en una gráfica que incluye en el eje X el tiempo ordenado cronológicamente y en el eje Y las variaciones de los datos unida por una linea y sin valores ausentes. Para esto se utiliza la función ts_plot, acompañada de los nombres de los ejes. Más información en: https://jesusbz.github.io/modelo-ventas.html

Figura 1. Ploteo de la serie de tiempo con los datos originales, para iniciar un análisis exploratorio de forma visual donde se podrá identificar si existe una tendencia en los datos de incremento o decremento y si presenta cambios en la variabilidad a travez del tiempo.

Componente estacional:

Permite identificar un patron estacional, mostrando la distibución de las muertes separada para cada año considerado, en cada mes. Esta figura es muy útil y posee un recurso interactivo interesante, ya que en la margen derecha, estan los colores que representan cada año y parandose en el que desen, pueden apagarlo o prenderlo, con el fin de verificar cada año por separado. Nuevamente si el cursor pasa por cualquier lines dibujada, se presentará el dato del mes y el valor, que en este caso es la cantidad de muertes totales correspondiente.

Figura 2. Imagen interactiva donde se muestra el componente estacional para cada año de forma independiente de la serie de tiempo

Componente estacional:

Se utiliza para identificar tendencias y patrones entre unidades del ciclo, en este caso son los meses los que se pueden manipular para ver las tendencias de cada uno o compararlos con otros de ellos, utilizando para esto la herramienta que aparece en la mano derecha, donde se pueden apagar o prender cada mes, según lo deseado.

Figura 3. Imagen interactiva donde se muestra el componente estacional para cada mes de forma independiente de la serie de tiempo

Componente estacional:

Se utiliza para identificar tendencias y patrones entre unidades del ciclo pero con box plot, que muestra la distribución de los datos en cada unidad del ciclo. Esta figura es muy importante para ver la distibución de los datos para cada mes en el rango de estudio. Nuevamente presenta la posibilidad de apagar o prender cada mes por separado y realizar comparaciones.

Figura 4. Imagen interactiva donde se muestra la distribución de datos en box plot por cada mes a lo largo de los 20 años considerados

Componentes totales

En la figura siguiente se muestra otro recurso visual pero en este caso tidimencional, donde esta representada la relación entre años, meses y valores en este caso de muertes totales. En esta figura aparece dos herramienta adicional, en la parte superior que permitirán rotar la imagen en cualquier eje, para ver con más claridad la distribución de los datos. Además, al desplazarse por la imagen se puede contrar en cualquier punto la informaci´pon de las tres variables consideradas aquí.

Si desea revisar más sobre el recurso visite: https://rpubs.com/ramkrisp/TSstudio

Figura 5. Imagen interactiva donde se muestra una gráfica tridimencional relacionando año, mes y cantidad de muertes totales en la serie de tiempo

Análisis visual de los datos buscando estacionariedad

Todos los análisis de series de tiempo deben cumplir con dos requisitos necesarios para poder confiar en el fin último de predecir valores a corto plazo, utilizando para esto datos históricos.

Los datos deben cumplir con la condición de estacionariedad, es decir que no tengan una tendencia de incremento o decremento, lo que es lo mismo decir: media cero a lo largo de su extensión
Cumplir con la condición de homocedasticidad, es decir que los datos se mantengan a lo largo del tiempo con la misma variabilidad.

Para esto hay varias formas de identificar estas condiciones. La primera es visual, donde se plotea la descomposición de los datos en tres componentes:

Serie observada = Tendencia + Efecto estacional + Residuos.

En la imagen aparecerá la figura de los datos completos (Observed), tendencia (Trend), que es el recurso para verificar si presenta incremento o decremento en la serie de tiempo, estacionalidad (Seasonal), donde se puede identificar si se presentan condiciones repetidas cada cierto tiempo) y los datos aleatorios (Random o residuos), que es el componente aleatorio de los datos debido a otras variables fortuitas. La función que se utiliza es: ts_plot(decompose(Datos.ts, type = ““)), se podría especificar en type, si se quiere additive (aditivo) i multiplicative (multiplicativo), o como en este caso los dos tipos al mismo tiempo (”both”), Si se detecta visualmente una heterocedasticidad de los datos, se tendría en cuenta el tipo multiplicative. Una de las transformaciones que podrían utilizarse cuando se detecte heterocedasticidad de los datos a lo largo del tiempo, es la utilización de logaritmo natural (log). https://www.youtube.com/watch?v=w-FAPJ9kLo4

ts_decompose(Datos.ts, type = "both")

Figura 6. Descomposición de la serie de tiempo para identificar los componentes de tendencia, estacionalidad y elementos aleatorios presentes.

El segundo recurso gráfico, hace referencia a lo que se conoce como figura para el análisis autregresivo, donde se establece una figura que presenta la correlación del valor en el momento t con el t-1 (un rezago) y de ahí hacia atrás, que puede dar información importante para tomar decisiones respecto a lo que se debería hacer con los datos antes del análisis. Para saber como se calculan la utocorrelación, sugiero que se vean erl siguiente tutorial: https://www.youtube.com/watch?v=ioRbWW-oEVI. Importante tener en cuenta que la gráfica auto regresiva hace referencia especialmente al componente de medias moviles, y cuantas se deberian tener en cuenta. Como se puede ver abajo, la disposición de la autoregresión primera, es típica de datos sin estacionariedad, es decir que los datos de t y t-1 se encuentan correlacionados altamente y van disminuyendo esa correlación, mientras más lejos esté históricamente, generando una reducción paulatina.

Para generar este recurso visual se utilizará el script: ts_cor(Datos.ts,lag.max = 60), el cual permitirá ver las autocorrelaciones y las autocorrelacions parciales con retrazo de 60 meses, y se podrán identificar posibles estacionalidades en la serie de tiempo.

La lineas punteadas horizontales y paralelas de color verde, estan mostrando la banda de confianza al 95%, donde se presentaria lo que se conoce como ruido blanco (que no presenta autocorrelación negativa o positiva)

ts_cor(Datos.ts,lag.max = 60)

Figura 7.Resultados de la autoregresión de los valores originales de la serie de tiempo. Se puede ver que se identifica la estacionalidad (barras rojas) cada 12 meses. Las linea verdes punteadas, delimitan una zona de no significancia de la correlación a un 95%.

Como se puede ver la figura de autocorrelación (parte superior), tiene ese descenso paulatino de las correlaciones típico de los datos con estacionariedad. Pero además identifica una estacionalidad cada 12 meses (barras rojas) y la figura de abajo, se muestra el modelo autoregresivo pero parcial, que identifica nuevamente la estacionalidad cada 12 meses.

Graficar por mes

Buscando estacionariedad en la estacionalidad, es decir que tambien debe presentarse media cero y homocedasticidad en cada unidad de tiempo, en este caso en cada mes a lo largo del año. El comando cycle determina la unidad de tiempo a la que pertenece cada observación de la serie, en este caso son 12 meses. Los resultados que se presentan es que existe para cada mes una tendencia creciente, pero parece existir homocedasticidad de los datos. https://rpubs.com/palominoM/series

par(mfrow=c(1,2), mar=c(4,4,4,1)+.1)
Datosmes<-monthplot(Datos.ts, col = "midnightblue",ylab="Datos originales", xlab= "Meses considerados en los 20 años")
Datosmes
Datos<-Datos.ts
boxplot(Datos~cycle(Datos))

Figura 8.Analisis de tendencia y homocedasticidad de los datos originales para cada unidad de tiempo considerada en este caso (mes), que permitirá tomar decisiones sobre la transformación de los datos para asegurar estacionariedad.

## Error : The fig.showtext code chunk option must be TRUE
## NULL
## Error : The fig.showtext code chunk option must be TRUE

Análisis estadísticos buscando estacionariedad

Esta es una forma más clara y contundente de identificar la estacionariedad de la serie de tiempo y es la utilización del test Dickey-Fuller, donde plantea:

Hipotesis nula (Ho): la serie no presenta estacionariedad p>0.05.

Hipotesis alterna (Ha): la serie presenta estacionariedad p<0.05

adf.test(Datos.ts, alternative="stationary", k=1)

## Warning in adf.test(Datos.ts, alternative = "stationary", k = 1): p-value
## smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Datos.ts
## Dickey-Fuller = -8.0784, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationary

El resultado encontrado, diria que la serie presentaria estacionariedad, ya que p es menor que 0.05. Sin embargo, al solicitar la cantidad de diferenciaciones que se requieren en los datos, se encuentra que es necesario un paso adicional

Transformación de datos:

Cuando se encuentra que no hay estacionariedad de los datos, es importante transformarlos para que cumplan esta condición, es así como se recomienda utilizar el siguiente script, que nos dirá el ndiffs: número de diferencias requeridas para que la serie total sea estacionaria. nsdiffs: determine the numero of diferencias estacionales se requerien para convertirlo en estacionario. El resultado de este análisis es que se requiere una diferenciación en los datos totales y no se necesita realizar ninguna diferenciación en los datos estacionales. En la siguiente figura se puede entender como se diferencian los datos con un retrazo y con dos, lo que en teoría convertirán datos estacionarios en datos no estacionarios

Cálculo que se realiza para encontrar la diferenciación para 1 y para 2 niveles.

A continuación se evaluan la cantidad de diferenciaciones necesarias que se requieren para convertir los datos totales con estacionariedad utilizando la función ndiffes y la cantidad de diferenciación requerida para eliminar la estacionariedd de la estacionalidad, utilizando la función nsdiffs, si es que existe.

ndiffs(Datos.ts)

## [1] 1

nsdiffs(Datos.ts)

## [1] 0

El resultado muestra que es necesario una diferenciación para los datos y que no se requiere diferenciación para eliminar la estacionariedad en la estacionalidad.

Transformación de datos:

La siguiente es una diferenciación para eliminar la tendencia, lo que convertirá a la serie de tiempo, para poder hacer análisis y predicciones. El siguiente script permite evidenciar la no estacionariedad de los datos ya transformados por la diferenciación de nivel 1. Si en el la linea 127 aparece que la diferenciación es 2 se utilizaria el siguiente comando: diff(x, lag = 2), lo que convierte a los datos reales en datos estacionales, que es lo que se busca.

seriedif1=diff(Datos.ts)
ts_plot(seriedif1)

Figura 9.Ploteo de los datos ya transformados teniendo en cuenta una diferenciación, donde se puede ver visualmente la posible estacionariedad de los datos despues de la transformación sugerida y la homocedasticidad.

Comparación datos originales (no estacionales) y de los datos con diferencias (estacionales)

En esta figura se puede apreciar las diferencias entre datos no estacionales y los datos estacionales. Cuando se revisa las dos imagenes de la columna izquierda las diferencias son evidentes. La gráfica de arriba muestra incremento a lo largo del tiempo, mientras que en la gráfica de abajo se ha eliminado ese incremento, con una media en 0, y la variabilidad no cambia visualmente, es decir que presenta homocedasticidad. En lo que respecta a las dos imagenes de la columna derecha, se presenta una comparación entre las graficas autoregresivas, donde en la parte superior presenta una típica de datos no estacionarios, es decir que los valores de autoregresión van disminuyendo paulatinamente, mientras se alejan de el valor más reciente, siendo el valor inicial relacionado con el mismo, mientras que la gráfica de abajo, presenta una autocorrelación de los valores estacionarios donde se elimina esta tendencia paulatina. Además, aparecen dos valores altamente relacionados: el primero de ellos negqativamente con un retrazo de un mes y el segundo positivamente con un retrazo de 12 meses.

par(mfrow=c(2,2), mar=c(4,4,4,1)+.1)
plot(Datos.ts, ylab= "Muertes totales")
acf(Datos.ts, main="Serie que no presenta estacionariedad")
plot(seriedif1)
acf(seriedif1, main="Serie que presenta estacionariedad")

Figura 10.Comparación de datos originales que no presentan estacionariedad y los datos transformados realizando una diferenciación, que si la presentan.Se muestran los datos de serie de tiempo a la izquierda y las gráficas autoregresivas correspondientes.

## Error : The fig.showtext code chunk option must be TRUE
## Error : The fig.showtext code chunk option must be TRUE
## Error : The fig.showtext code chunk option must be TRUE
## Error : The fig.showtext code chunk option must be TRUE

Graficar por mes

https://rpubs.com/palominoM/series Este recurso es muy importante porque puede determinar si dentro de los meses tambien se presenta estacionariedad, es decir que no se presente una tendencia a incrementar o decrementar los valores (media=0) y si presenta homocedasticidad que es lo que se quiere. La gáfica de la izquierda muestra que no hay un incremento o decremento (no tiene tendencia), y la figura de la derecha muestra la homocedasticidad para cada mes.

par(mfrow=c(1,2), mar=c(4,4,4,1)+.1)
Estames<-monthplot(seriedif1, col = "midnightblue")
Estames
Muertesmes<-seriedif1
boxplot(Muertesmes~cycle(Muertesmes))

Figura 11.La gráfica a la izquierda muestra la evidencia que no hay una tendencia de los datos ya transformados de serie de tiempo para cada mes por separado y la gráfica de la derecha, muestra la homocedasticidad de los datos tambien para cada mes.

## Error : The fig.showtext code chunk option must be TRUE
## NULL
## Error : The fig.showtext code chunk option must be TRUE

Autocorrelacion parcial

Mide la correlaciíon existente entre una variable en distintos periodos de tiempo, pero eliminando el efecto de periodos intermedios entre el rango escogido. Si se desea más información del calculo, puede visitar este tutorial: https://www.youtube.com/watch?v=GLbV_VEkZVA Importante tener en cuenta que la gráfica de autoregresivos parciales hace referencia especialmente al componente autoregresivo y cuantos autoregresivos se deberian tener en cuenta. Para identificar en el modelo cuantos autoregresivos se debe considerar con el modelo original https://www.youtube.com/watch?v=PBcdIK5oKtQ&t=8s 14’ Cuando se realiza una autoregreción parcial, nos ayuda a entender algunas cosas de los datos trabajados. Por ejemplo cuando la primera relación sea significativa y abruptamente se reduzca a valores no significativos, estariamos posiblemente frenta a valores no estacionarios, que requieren por ejemplo una diferenciación. , puede determinarse los rezagos que hay que tener en cuenta en el modelado de los datos.

ts_cor(seriedif1, lag.max =36)

Figura 12.Comparación de la autoregresión parcial de los datos originales a la izquierda y los datos transformados a la derecha.

Figura de correlaciones con retrazos

En esta figura se muestran las autocorrelaciones que se presentan, con retrazo de 1, 6 y 12 meses. En estos tres retrazos se pueden ver que el primero tiene una correlación negativa, el de 6 meses continua siendo negativa, péro el de 12 meses es claramente positiva

ts_lags(seriedif1, lags = c(1,6,12))

Figura 13.Comparación y relación entre autocorrelaciones en los rezagos 1, 6 y 12 meses, para la serie de tiempo.

Modelos de ajuste:

El modelo que utilizaremos es el conocido como ARIMA, que se descompone en los siguientes elementos

AR (Componente autoregresivo): hace referencia a la relación que se establece entre los datos y ellos mismos, que permitirá predecir lo que pasa hoy con respecto a lo que paso ayer o considerando un mayor rezago.

I (Componente de diferenciación): que puede ser de diferentes ordenes (1, 2, 3, etc.) y que se utilizan para generar estacionariedad en la serie de tiempo, que es una condición para trabajar en este tema. Ver: https://www.youtube.com/watch?v=FeohLw3dxJE (No se trabaja con los valores reales sino con los cambios con respecto al t-1

MA (Medias moviles): tiene que ver con los errores (ruido blanco)

S (Estacionalidad): fenómeno que se repite

Parámetros del modelo: (p,d,q)(ps,ds,qs)[S] P: El orden del AR

D (diferencias): Este modelo permitirá identificar los elementos que se hizo referencia arriba para ajustarse a los datos de mejor forma.

En los resultados se puede identificar los parametros del modelo, que son los valores que se encuentran entre parentesis frente al termino ARIMA. El primer parentesis (parte regular del modelo) tiene valores (1,0,1), que estan representando en su orden el componente autoregresivo, componente de diferencición, el cual aparece 0 ya que estamos trabajando con la serie de tiempo diferenciada y componente de medias moviles, que corresponden a los primeros coeficientes (ar1 y ma1). En cuanto se refiere al parentesis segundo (0,0,2)(parte estacional), solo aparece el componente de medias moviles y corresponde a los otros dos coeficientes, en el resultado despues de correr el modelo. Importante decir que si se tienen varios modelos, la decisión para escoger el mejor ajuste, se utilizaría el AIC que genera el modelo o el valor log likelihood, siempre favoreciendo el valor menor.

ModeloARIMA<- auto.arima(seriedif1)
ModeloARIMA

## Series: seriedif1 
## ARIMA(1,0,1)(0,0,2)[12] with non-zero mean 
## 
## Coefficients:
##          ar1      ma1    sma1    sma2     mean
##       0.3402  -0.9413  0.4315  0.2005  23.6392
## s.e.  0.0908   0.0493  0.0708  0.0588   8.1166
## 
## sigma^2 = 646744:  log likelihood = -1937.33
## AIC=3886.66   AICc=3887.02   BIC=3907.52

Supuestos del modelo:

Relacionado con los residuos que deben cumplir con algunas caracterìsticas importantes. Estos residuos parecen ruido blanco, por lo que su serie presenta estacionariedad y modelada correctamente. No existe autocorrelación de los residuos.

tsdiag(ModeloARIMA)

Figura 13. Presentación de los supuestos de los residuales del modelo auto.arima, donde se puedem ver que cumple todos ellos.

## Error : The fig.showtext code chunk option must be TRUE
## Error : The fig.showtext code chunk option must be TRUE
## Error : The fig.showtext code chunk option must be TRUE

Diagnosis del modelo Metodo 1

En este diagnostico tenemos:

Ho: Los residuos del modelo presentan autocorrelación p<0.05

Ha: Los residuos del modelo presentan no presentan autocorrelación p>0.05

En este caso los residuos no presentan autocorrelación, es decir que el modelo es correcto, y cumple con los supuestos requeridos.

Box.test(ModeloARIMA$residuals, lag=10, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  ModeloARIMA$residuals
## X-squared = 5.1746, df = 10, p-value = 0.8792

Diagnosis del modelo Metodo 2

La gráfica superior no muestra una tendencia y varia al rededor del 0 de forma aleatoria, que es lo que se espera en estos casos. En la gráfica inferior izquierda, se encuentra una grafica autoregresiva, con las lineas azules punteadas paralelas, que determinan una zona donde se esperaria que los valores de los residuos no sobrepasen. En este caso solo hay uno que sobrepasa una de las lineas puneteadas, pero el valor de correlación es muy bajo. Y finalmente en la gáfica inferior derecha, se presenta la distribución normal de los residuos, que a pesar de tener una cola larga hacia la derecha, tambien se acerca a la distribución normal. https://towardsdatascience.com/f-forecasting-5d23341462eb

autoplot(ModeloARIMA,seasonal = TRUE)

Figura 14. Presentación del método 2 de diagnostico, donde se puede ver el cumplimiento ellos.

check_res(ModeloARIMA)

## Error : The fig.showtext code chunk option must be TRUE

Figura 14. Presentación del método 2 de diagnostico, donde se puede ver el cumplimiento ellos.

Comparacion de los datos observados con los esperados

En la siguiente figura se puede revisar el ajuste del modelo de forma visual

par(mfrow=c(1,1))
plot(ModeloARIMA$x, col="red")
lines(fitted(ModeloARIMA), col="blue")

Figura 15.Comparación visula del ajuste del modelo (Rojo) y los datos reales tranformados (azul).

## Error : The fig.showtext code chunk option must be TRUE

Predicción (Pronóstico)

La siguiente figura presenta la predicciòn con el modelo autoarima que se genero despues que se transformaron los datos segun lo necesario para asegurar estacionariedad y homocedasticidad. Los datos que aparecen a continuación es el pronostico que se consigue con el modelo autoarima, teniendo en cuenta que la predición se hace a 12 meses adelante, es decir que corresponde a 1999 de enero a diciembre y presenta los intervalos de confianza. Para comparar modelos revisar minuto 20 https://www.youtube.com/watch?v=a5QQp9peaZ4

forecast(ModeloARIMA, h=12)

##          Point Forecast      Lo 80     Hi 80     Lo 95    Hi 95
## Jan 1999     -444.10019 -1474.7294  586.5290 -2020.312 1132.111
## Feb 1999     -812.99812 -2015.5110  389.5148 -2652.083 1026.087
## Mar 1999      377.26846  -843.5732 1598.1101 -1489.848 2244.385
## Apr 1999     -354.66909 -1577.6140  868.2758 -2225.002 1515.664
## May 1999      502.61134  -720.5767 1725.7994 -1368.093 2373.316
## Jun 1999     -119.06020 -1342.2764 1104.1559 -1989.808 1751.687
## Jul 1999      164.19170 -1059.0277 1387.4111 -1706.561 2034.944
## Aug 1999      148.25476 -1074.9650 1371.4745 -1722.498 2019.008
## Sep 1999      -53.37802 -1276.5978 1169.8418 -1924.131 1817.375
## Oct 1999       35.73000 -1187.4898 1258.9498 -1835.023 1906.483
## Nov 1999       53.80178 -1169.4180 1277.0216 -1816.952 1924.555
## Dec 1999      320.98009  -902.2397 1544.1999 -1549.773 2191.733

Comparación visual de los pronosticos

Se genera el modelo autoarima con los datos originales

Modeloarima<-auto.arima(Datos.ts)
Modeloarima

## Series: Datos.ts 
## ARIMA(1,1,1)(0,0,2)[12] with drift 
## 
## Coefficients:
##          ar1      ma1    sma1    sma2    drift
##       0.3402  -0.9413  0.4315  0.2005  23.6392
## s.e.  0.0908   0.0493  0.0708  0.0588   8.1166
## 
## sigma^2 = 646745:  log likelihood = -1937.33
## AIC=3886.66   AICc=3887.02   BIC=3907.52

En este caso lo que se pretende es mostrar el modelo autoarima que se consiguio con los valores históricos de los datos reales y los de los datos diferenciados y el pronóstrico para un año adelante.

par(mfrow=c(2,1), mar=c(4,4,4,1)+.1)
autoplot(forecast::forecast(Modeloarima,h=12))

Figura 15.Comparación visula del ajuste del modelo (Rojo) y los datos reales tranformados (azul).

autoplot(forecast::forecast(ModeloARIMA,h=12))

Figura 15.Comparación visula del ajuste del modelo (Rojo) y los datos reales tranformados (azul).

## Error : The fig.showtext code chunk option must be TRUE
## Error : The fig.showtext code chunk option must be TRUE

Modelo autoarina para los datos originales

Modelo01<-auto.arima(Datos.ts)
Modelo01

## Series: Datos.ts 
## ARIMA(1,1,1)(0,0,2)[12] with drift 
## 
## Coefficients:
##          ar1      ma1    sma1    sma2    drift
##       0.3402  -0.9413  0.4315  0.2005  23.6392
## s.e.  0.0908   0.0493  0.0708  0.0588   8.1166
## 
## sigma^2 = 646745:  log likelihood = -1937.33
## AIC=3886.66   AICc=3887.02   BIC=3907.52

Evalución visual del modelo

Se grafican los valores reales con los valores justados con el modelo auto.arima

plot(Modelo01$x, col="red",xlab="Tiempo (Año y meses)",ylab="Muertes totales")
lines(fitted(Modelo01), col="blue")

Figura 15.Comparación visula del ajuste del modelo (Rojo) y los datos reales tranformados (azul).

## Error : The fig.showtext code chunk option must be TRUE

Evaluación de la eficacia del modelo

Se podría establecer un modelo de correlación entre los datos reales del rango entre 1979 a 1998 y los datos ajustados del modelo, para tener un valor de correlación (R2), que podría ser una medida de eficacia. En este caso el valor fue de 0.78, es decir que los valores de reales son pronosticados por los valores ajustados en un 78%.

Eficaciamodelo<-lm(Modelo01$x~fitted(Modelo01))
summary(Eficaciamodelo)

## 
## Call:
## lm(formula = Modelo01$x ~ fitted(Modelo01))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2441.7  -387.2   -27.1   313.6  6307.0 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      854.66988  416.39062   2.053   0.0412 *  
## fitted(Modelo01)   0.93311    0.03239  28.809   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 790.4 on 238 degrees of freedom
## Multiple R-squared:  0.7771, Adjusted R-squared:  0.7762 
## F-statistic: 829.9 on 1 and 238 DF,  p-value: < 2.2e-16

Otra forma de evaluar el modelo encontrado sería comparar los valores predichos por el modelo y los valores reales del siguiente año, sea porque los tenemos en una base de datos o sea porque esperamos la finalizaciòn de ese nuevo año y lo comparamos con lo predicho hace una año. Esta sería una medida muy objetiva de la bondad de ajuste del modelo seleccionado. Es importante entender que cuando se hace un pronóstico, se debería restringir a tiempos muy cortos, por ejemplo a 12 meses en este caso. Es así que si pudieramos contar con los datos reales de ese año, en este caso 1999, ya que nuestro modelo se base en datos entre 1979 y 1998., se podrìa establecer una evaluación y generar ajustes al modelo, adicionando esa informaciòn.

Se obtienen los valores reales del año pronosticado por el modelo, es decir los muertos totales de 1999 para cada mes

Colombia99 <- read_excel("Col99.xlsx")

## New names:
## • `` -> `...1`

head(Colombia99)

## # A tibble: 6 × 3
##   ...1  MES   conteo
##   <chr> <chr>  <dbl>
## 1 1     Ene    16015
## 2 2     Feb    13287
## 3 3     Mar    14902
## 4 4     Abr    14191
## 5 5     May    15570
## 6 6     Jun    15382

Aquì se organiza la informaciòn predica de los doce meses de 1999, utilizando el modelo ARIMA, utilizando para esto el modelo con los datos originales

X<-forecast(Modelo01, h=12)
Predichos<-as.data.frame(X)
Predichos<-Predichos[,-c(2:5)]
Predichos<-as.data.frame(Predichos)
MES<- c("Ene", "Feb", "Mar", "Abr", "May", "Jun", "Jul", "Ago", "Sep", "Oct", "Nov", "Dic")
MES<-as.data.frame(MES)
Predicho99<-cbind(MES,Predichos)

Ajuste<-merge(Predicho99,Colombia99, by="MES")

ggplot(Ajuste, aes(x=conteo, y=Predichos)) + 
  geom_point() + theme_light()+ geom_smooth(method='lm', formula=y~x, se=FALSE, col='dodgerblue1')

## Error : The fig.showtext code chunk option must be TRUE

mod1 <- lm(Predichos ~ conteo, data=Ajuste)
summary(mod1)

## 
## Call:
## lm(formula = Predichos ~ conteo, data = Ajuste)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -344.24 -106.32  -59.85  117.23  405.53 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.146e+04  1.139e+03  10.056 1.51e-06 ***
## conteo      2.312e-01  7.430e-02   3.112    0.011 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 234.7 on 10 degrees of freedom
## Multiple R-squared:  0.492,  Adjusted R-squared:  0.4412 
## F-statistic: 9.684 on 1 and 10 DF,  p-value: 0.01103

Algoritmo machine learning para series de tiempo (Entrenar 80% y testeo 20%)

Libreria necesaria primer método

library(xgboost)

## 
## Attaching package: 'xgboost'

## The following object is masked from 'package:dplyr':
## 
##     slice

library(tidymodels)

## ── Attaching packages ────────────────────────────────────── tidymodels 1.0.0 ──

## ✔ broom        1.0.1     ✔ rsample      1.1.0
## ✔ dials        1.0.0     ✔ tune         1.0.1
## ✔ infer        1.0.3     ✔ workflows    1.1.0
## ✔ modeldata    1.0.1     ✔ workflowsets 1.0.0
## ✔ parsnip      1.0.2     ✔ yardstick    1.1.0
## ✔ recipes      1.0.2

## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ yardstick::accuracy() masks forecast::accuracy()
## ✖ scales::discard()     masks purrr::discard()
## ✖ dplyr::filter()       masks mice::filter(), stats::filter()
## ✖ recipes::fixed()      masks stringr::fixed()
## ✖ dplyr::lag()          masks stats::lag()
## ✖ recipes::prepare()    masks VIM::prepare()
## ✖ xgboost::slice()      masks dplyr::slice()
## ✖ yardstick::spec()     masks readr::spec()
## ✖ recipes::step()       masks stats::step()
## • Dig deeper into tidy modeling with R at https://www.tmwr.org

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(timetk)

Convertir la serie de tiempo en dataframe y agregarle una columna de fecha

https://www.youtube.com/watch?v=LZu9LHQFNZw&t=1461s db<-dataset %>% dplyr::mutate(fecha_corte::tsibble::yearmonth(fecha_corte))

DatosT<-as.data.frame(Datos.ts)
f <- as.Date("1979-01-01")
DatosT$Fecha<-seq(from=f, by="months", length.out=240)
require(reshape)

## Loading required package: reshape

## 
## Attaching package: 'reshape'

## The following object is masked from 'package:lubridate':
## 
##     stamp

## The following object is masked from 'package:dplyr':
## 
##     rename

## The following objects are masked from 'package:tidyr':
## 
##     expand, smiths

DatosT = rename(DatosT, c(x="Valores"))

Visualización interactiva

DatosT %>% 
  plot_time_series(Fecha,DatosT$Valores,.interactive = TRUE)

## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.

Generar los datos de entrenamiento y testeo del modelo

Importante tener en cuenta que cuando se generan los datos para entrenamiento y testeo, no se puede realizar aleatoriamente, sino deben separarse dos grupos sin perder la cronología de los datos, es decir que el 80% seran los datos primeros y el 20% deberian ser los datos finales

splits= initial_time_split(DatosT, prop=0.9)
glimpse(training(splits))

## Rows: 216
## Columns: 2
## $ Valores <dbl> 9582, 8467, 9202, 9197, 9484, 9596, 9560, 9049, 8424, 8988, 87…
## $ Fecha   <date> 1979-01-01, 1979-02-01, 1979-03-01, 1979-04-01, 1979-05-01, 1…

plot.ts(training(splits))

## Error : The fig.showtext code chunk option must be TRUE

## Error : The fig.showtext code chunk option must be TRUE

glimpse(testing(splits))

## Rows: 24
## Columns: 2
## $ Valores <dbl> 15042, 12800, 14204, 13892, 14699, 14403, 14706, 14430, 13536,…
## $ Fecha   <date> 1997-01-01, 1997-02-01, 1997-03-01, 1997-04-01, 1997-05-01, 1…

plot.ts(testing(splits))

## Error : The fig.showtext code chunk option must be TRUE

## Error : The fig.showtext code chunk option must be TRUE

Entrenar el modelo (auto-arima)

Modelo1Nboots<-arima_reg() %>% 
  set_engine(engine = "auto_arima") %>% 
  fit(Valores~Fecha,data=training(splits))

## frequency = 12 observations per 1 year

Modelo1Nboots

## parsnip model object
## 
## Series: outcome 
## ARIMA(1,1,1)(0,0,2)[12] with drift 
## 
## Coefficients:
##          ar1      ma1    sma1    sma2    drift
##       0.3319  -0.9418  0.4275  0.1862  25.3368
## s.e.  0.1022   0.0591  0.0754  0.0631   8.6255
## 
## sigma^2 = 682671:  log likelihood = -1748.48
## AIC=3508.96   AICc=3509.36   BIC=3529.18

Calibracion del modelo

Calibracion<-Modelo1Nboots %>% 
  modeltime_calibrate(new_data = testing(splits))

## Converting to Modeltime Table.

Calibracion$.calibration_data

## [[1]]
## # A tibble: 24 × 4
##    Fecha      .actual .prediction .residuals
##    <date>       <dbl>       <dbl>      <dbl>
##  1 1997-01-01   15042      14884.      158. 
##  2 1997-02-01   12800      14338.    -1538. 
##  3 1997-03-01   14204      14748.     -544. 
##  4 1997-04-01   13892      14599.     -707. 
##  5 1997-05-01   14699      14757.      -58.0
##  6 1997-06-01   14403      14736.     -333. 
##  7 1997-07-01   14706      14869.     -163. 
##  8 1997-08-01   14430      14820.     -390. 
##  9 1997-09-01   13536      14904.    -1368. 
## 10 1997-10-01   13604      15207.    -1603. 
## # … with 14 more rows

Ploteo real y predicho

Calibracion %>% 
  modeltime_forecast(new_data = testing(splits),actual_data = DatosT) %>% 
  plot_modeltime_forecast(.legend_max_width = 25, .interactive = T)

Evaluacion del o los modelo(s)

Modelos adicionales par agregar: https://www.youtube.com/watch?v=iC4Z1ldNRxY&t=3967s En 1h 8’ se adiciona posibilidades de manejar manualmente el orden=c (0,1,1) y estacionalidad seasonal=c(0,0,1)

Donde se mira estos indices: https://www.youtube.com/watch?v=5SJrPLRJJ8Q&t=2322s 38’ MAE: Error absoluto medio

MAPE: Porcentaje de error absoluto medio

MASE: Escala de error absoluto medio

SMAPE: Porcentaje de error abosouto medio

RMSE: Root mean squared error

RSQ: R cuadrado

Calibracion %>% 
  modeltime_accuracy() %>% 
  table_modeltime_accuracy(.interactive = T)

Libreria necesaria segundo método

Revisar: https://github.com/RamiKrispin/TSstudio

library(TSstudio)
library(forecast)

Modelo de pronostico de entrenamiento

Dtos_s <- ts_split(ts.obj = Datos.ts, sample.out = 12)
train <- Dtos_s$train
test <- Dtos_s$test

Predicción con autoarima

md <- auto.arima(train)
fc <- forecast(md, h = 12)

Ploteo datos reales y ajustados y la predicción

test_forecast(actual = Datos.ts, forecast.obj = fc, test = test)

Ploteo de predicción

plot_forecast(fc)

Comparación de diferentes métodos

methods <- list(ets1 = list(method = "ets",
                            method_arg = list(opt.crit = "lik"),
                            notes = "ETS model with opt.crit = lik"),
                ets2 = list(method = "ets",
                            method_arg = list(opt.crit = "amse"),
                            notes = "ETS model with opt.crit = amse"),
                arima1 = list(method = "arima",
                              method_arg = list(order = c(1,1,0)),
                              notes = "ARIMA(1,1,1)"),
                arima2 = list(method = "arima",
                              method_arg = list(order = c(1,1,1),
                                                seasonal = list(order = c(0,0,2))),
                              notes = "SARIMA(1,1,1)(0,0,2)"),
                hw = list(method = "HoltWinters",
                          method_arg = NULL,
                          notes = "HoltWinters Model"),
                tslm = list(method = "tslm",
                            method_arg = list(formula = input ~ trend + season),
                            notes = "tslm model with trend and seasonal components"))

Entrenando los modelos

md <- train_model(input = Datos.ts,
                  methods = methods,
                  train_method = list(partitions = 6, 
                                      sample.out = 12, 
                                      space = 3),
                  horizon = 12,
                  error = "MAPE")

## Warning in (function (x, alpha = NULL, beta = NULL, gamma = NULL, seasonal =
## c("additive", : optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH

## # A tibble: 6 × 7
##   model_id model       notes                     avg_m…¹ avg_r…² avg_c…³ avg_c…⁴
##   <chr>    <chr>       <chr>                       <dbl>   <dbl>   <dbl>   <dbl>
## 1 ets1     ets         ETS model with opt.crit …  0.0232    432.   1           1
## 2 hw       HoltWinters HoltWinters Model          0.0283    494.   1           1
## 3 ets2     ets         ETS model with opt.crit …  0.0293    498.   0.958       1
## 4 arima2   arima       SARIMA(1,1,1)(0,0,2)       0.0340    612.   0.944       1
## 5 arima1   arima       ARIMA(1,1,1)               0.0451    792.   0.972       1
## 6 tslm     tslm        tslm model with trend an…  0.0491    768.   0.903       1
## # … with abbreviated variable names ¹avg_mape, ²avg_rmse, ³`avg_coverage_80%`,
## #   ⁴`avg_coverage_95%`

Ploteo del rendimiento de los diferentes modelos en las particiones de prueba

plot_model(md)

Series de tiempo

Oscar R Sánchez

25/11/2022

Introducción

Descargar los datos:

Descarga de base de datos:

Revisión de base de datos

Datos faltantes

Librerias necesarias:

Ejemplo de imputación

Librerias necesarias:

Conversión a serie de tiempo:

Análisis exploratorio :

Componente estacionario

Componente estacional:

Componente estacional:

Componente estacional:

Componentes totales

Análisis visual de los datos buscando estacionariedad

Graficar por mes

Análisis estadísticos buscando estacionariedad

Transformación de datos:

Transformación de datos:

Comparación datos originales (no estacionales) y de los datos con diferencias (estacionales)

Graficar por mes

Autocorrelacion parcial

Figura de correlaciones con retrazos

Modelos de ajuste:

Supuestos del modelo:

Diagnosis del modelo Metodo 1

Diagnosis del modelo Metodo 2

Comparacion de los datos observados con los esperados

Predicción (Pronóstico)

Comparación visual de los pronosticos

Evalución visual del modelo

Evaluación de la eficacia del modelo

Algoritmo machine learning para series de tiempo (Entrenar 80% y testeo 20%)

Libreria necesaria primer método

Convertir la serie de tiempo en dataframe y agregarle una columna de fecha

Visualización interactiva

Generar los datos de entrenamiento y testeo del modelo

Entrenar el modelo (auto-arima)

Calibracion del modelo

Ploteo real y predicho

Evaluacion del o los modelo(s)

Libreria necesaria segundo método

Modelo de pronostico de entrenamiento

Predicción con autoarima

Ploteo datos reales y ajustados y la predicción

Ploteo de predicción

Comparación de diferentes métodos

Entrenando los modelos

Ploteo del rendimiento de los diferentes modelos en las particiones de prueba